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Abstract: The truly meshless local Petrov-Galerkin 

(MLPG) method holds a great promise in solving boundary 
value problems, using a local symmetric weak form as a natu- 
ral approach. In the present paper, in the context of MLPG and 
the meshless interpolation of a moving least squares (MLS) 
type, a method which uses primary and secondary nodes in the 
domain and on the global boundary is introduced, in order to 
improve the accuracy of solution. The secondary nodes can be 
placed at any location where one needs to obtain a better res- 
olution, The sub-domains for the shape functions in the MLS 
approximation are defined only from the primary nodes, and 
the secondary nodes use the same sub-domains. The shape 
functions based on the MLS approximation, in an integration 
domain, have a single type of a rational function, which re- 
duces the difficulty of numerical integration to evaluate the 
weak form. The present method is very useful in an adap- 
tive calculation, because the secondary nodes can be easily 
added and/or moved without an additional mesh. The essential 
boundary conditions can be imposed exactly, and non-convex 
boundaries can be treated without special techniques. Several 
numerical examples are presented to illustrate the performance 
of the present method. 

keyword: meshiess method, MLPG method, local symmet- 

ric weak form, MLS, primary node, secondary node. 

1 Introduction 

Meshless methods are attractive in adaptive error-control in 
computations to solve boundary value problems, by adding or 
removing nodes without the burdensome task of remeshing, 
each time. Several meshless methods have been proposed, 
each with certain advantages and disadvantages. These in- 
clude: the smooth particle hydrodynamics (SPH) [Gingold and 
Monaghan (1977)], the element free Galerkin (EFG) method 
[Nayroles, Touzot, and Villon (1992)], the reproducing kernel 
particle method (RKPM) [Liu, Jun, and Zhang (1995)], the 
hp-cloud method [Duarte and Oden (1996)], the finite point 
method [Onate, Idelsohn, Zienkiewicz, and Taylor (1996)], 
the partition of unity [Babuska and Melenk (1997)], the lo- 
cal boundary integral equation (LBIE) method [Zhu, Zhang, 
and Atluri (1998a,b)], and the meshless local Petrov-Galerkin 
(MLPG) method [Atluri and Zhu (1998a,b)]. In these meth- 
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ods, the construction of a trial approximation, which does 
not rely on element connectivity, is a significant development. 
However, most meshless methods, except the LBIE/MLPG 
methods, are not truly meshless approaches, since these meth- 
ods require background meshes for numerical integration of 
the weak form. In these methods, an additional cost is asso- 
ciated with the construction of a background mesh, if nodes 
are added or deleted in a domain. The MLPG/LBIE methods, 
however, are more natural approaches, because these meth- 
ods use a local weak form, and use numerical integration over 
sub-domains, which can be of arbitrary shapes such as circles, 
ellipses, rectangulars and parallelopipeds in a 2-dimensional 
geometry. 

In spile of the novel concepts embodied in the MLPG method, 
difficulties in the numerical integration for evaluation of the 
weak form still persist, as reported by Atluri, Cbo, and Kim 
(1999a) and Atluri, Kim, and Cho (1999b). This is due to the 
complexity of the non-element interpolation functions, which 
result from the moving least squares (MLS), the partition of 
unity, and the hp-cloud methods. In addition, circular sub- 
domains make the numerical integration difficult, because the 
intersections between such sub-domains result in highly com- 
plex functions in the integration domain. As a result, a large 
number of integration points may be required to obtain ac- 
curate solutions [Atluri, Cho, and Kim (1999a) and Atluri, 
Kim, and Cho (1999b)]. In this paper, we present a viable 
method, based on the MLPG, that use secondary nodes to ob- 
tain an improvement in the accuracy of solution, without an 
additional mesh. The sub-domains for the MLS shape func- 
tions are generated only from the primary nodes, and the sec- 
ondary nodes use the same sub-domains. The secondary nodes 
in the domain, and on the global boundary, do not necessi- 
tate the creation of new sub-domains, and the shape functions 
for the secondary nodes can be easily defined on the origi- 
nal sub-domains, using the MLS approximation. Numerical 
integration is carried out over a polygonal cell, which is the 
intersection of the sub-domains constructed only from the pri- 
mary nodes. As a consequence of the alignment of the bound- 
aries of sub-domains and integration domains, the shape func- 
tions have a single type of a rational function in a domain of 
integration. This greatly alleviates the difficulty in the nu- 
merical integration of the weak form in the MLPG method. 
The essential boundary conditions can be imposed exactly, 
and the non-convex boundaries can be treated without using 
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Figure 1 : A schematic representation of the. sub-domain Q*, 
with node / as its center, and 9Q' as its boundary. The globa] 
domain is Q, with a global boundary T, where displacements 
are prescribed at T M) and tractions are prescribed at T r . 


any special techniques. The primary advantage of the present 
MLPG method is that the secondary nodes can be easily added 
and/or removed, without the burdensome task of constructing 
a new mesh, because the secondary nodes use the same sub- 
domains defined from the original primary nodes. The use of 
the secondary nodes to improve the computational solutions to 
a boundary value problem is independent of the mesh derived 
from the primary nodes. The present MLPG method offers 
a very useful tool for an adaptive calculation, by controlling 
errors in the computed results. 

There have been several efforts to develop ways to improve 
the accuracy of a numerical solution, using a coarse primary 
mesh. Oden, Duarte, and Zienkiewicz (1998) introduced a 
new hp-finite element method, by using a combination of the 
conventional FEM and the partition of unity, to achieve a dif- 
ferent order of basis for each node. In this method, how- 
ever, a new global mesh is needed to add nodes for refine- 
ments, and a careful choice of the basis functions has to be 
made to prevent their linear dependence. The so-called gen- 
eralized finite element (GFEM) [Strouboulis, Babuska, and 
Copps (1998)] uses special functions from known analytical 
solutions in order to improve the FEM solution, in a way that 
is similar to the conventional hybrid FEM [Atluri, Gallagher, 
and Zienkiewicz (1983)]. Liu, Uras, and Chen (1998) used a 
coupling of the RKPM and the FEM to achieve an adaptive 
calculation by adding nodes. Although this coupling of the 
RKPM and the FEM has features that are similar to the present 
MLPG method, the basic approaches of two methods are quite 
different. The enrichment using the RKPM is not based on a 
local weak-formulation, and the coupling of the RKPM and 
the FEM may not give rise to consistent solutions, because of 
the difficulty of numerical integration over an integration do- 
main when the boundaries of sub-domains and the integration 
domains are not aligned with each other. 


Figure 2 ; A schematic illustrating various shapes for sub- 
domains, and the region bounding all the nodes in Q which 
may influence the interpolation at a generic x in a meshless 
approximation. 

The outline of this paper is as follows. In section 2, the local 
symmetric weak form is explained, as a key concept in the 
MLPG method. In section 3, we review the characteristics of 
the MLS approximation and of the numerical integration. In 
this part, we emphasis the difficulty of numerical integration 
to evaluate the weak form in meshless methods. In section 4, 
the concept of ’’primary” and ’’secondary” nodes is introduced. 
To construct proper shape functions for the primary and the 
secondary nodes in the domain and on the global boundary, 
weight functions in the MLS approximation should be defined 
appropriately over sub-domains, in order to preserve a single 
type of a rational function in each integration domain. Towards 
this end, we present a method to construct the weight functions 
in the MLS approximation, for the primary and the secondary 
nodes. Numerical examples in linear elasticity are presented 
in section 5. Finally, the conclusions are given in section 6. 

2 The meshless local Petrov-Galerkin (MLPG) formula- 
tion 

The equilibrium equations of linear elasticity, in a global do- 
main Q bounded by T, are given by 

Oijj+bi = 0 in Cl (1) 

where G,y is the stress tensor, b x are the body forces, () j de- 
notes 0()/3xy, and a summation over a repeated index is im- 
plied. The boundary conditions are assumed to be 

u i ~ Tij at T u (2) 

o (J n j — 7, ai F, (3) 

where T u and r, are the global boundaries with prescribed 
displacements and tractions, respectively. In a conventional 
Galerkin finite element formulation, the global weak form is 
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used to solve the boundary value problem numerically. How- 
ever, the MLPG method starts from a weak form over a lo- 
cal sub-domain, or a patch, inside the global domain Q as 
shown in Fig. 1. Let {ft*] be a system of overlapping patches 
which covers the global domain ft, where /(= 1,2,*’- , TV) in- 
dicates a node, and N is the total number of nodes. We implic- 
itly introduce the concept of "nodes” with "local domains”. 
The sub-domain ft{ is thus called the sub-domain of node I. 
The sub-domain ft{ can be a circle, a rectangle, or an ellipse 
in two dimensions (or a sphere, a parallelopiped, or an ellip- 
soid in three dimensions) in the MLPG formulation, but it can 
be extended to any kinds of geometry as shown in Fig. 2. 

A generalized local weak form of the equilibrium equation is 
written as 


stiffness matrix, involving all the J nodes, whose sub-domains 
ft^ intersect with ft*, such that the integrand in Eq. 6 is non- 
zero. 

To obtain the discrete equations from the MLPG formulation 
in Eq. 6, based on meshless interpolations, which are ex- 
plained in the next section, the following interpolations are 
used. The global forms of interpolations for the trial and the 
test functions, respectively, can be written as 

= X (7) 

j= i 

v?(x) = (8) 

/=! 


J ( Oij J + b, ) Vi dQ = 0 (4) 

where v, is the test function. Using the divergence theorem in 
Eq. 4, we obtain the following local weak-form: 

f OijnjVjdT — [ (GijVij + bjVj)dCl = 0 (5) 

where is the outward unit normal to the boundary 3ft{, 

In the MLPG method, the trial and lest functions can be cho- 
sen from different spaces. In particular, the test functions need 
not vanish on the boundary where the essential boundary con- 
ditions are specified. Also, in order to simplify Eq. 5, we 
deliberately select the test functions v, such that they vanish 
over except when 3ft* intersects with the global bound- 
ary T. This can be easily accomplished in the MLPG method 
by using the nodal-tesl-shape function whose value at the local 
boundary 3ft{ is zero, as long as 3ft{ does not intersect with 
T. Using these test functions in Eq. 5, we obtain the following 
Local Symmetric Weak Form (LSWF): 


where § J {x) and \j / 7 (x) are the nodal shape functions for the 
trial and the test functions centered at nodes J and /, respec- 
tively. In general, in meshless interpolations, v{ and u J ( are 
fictitious nodal values. Substitution of Eq. 6 into the MLPG 
formulation in Eq. 5 leads to following discretized system of 
linear equation: 


y [ (Bl) T DB J it J dn= [ v't dr+ [ v'bdn 
jtj nj M, Jal 

where, in two-dimensional space, 
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[ GijVijdO.- [ //v/</r+ [ biVjdCl (6) 

M Jn's 

where r,- = Oijtij, and T* r is the intersection of T, and the 
boundary 3ft* of ft*. The MLPG method based on the lo- 
cal formulation in Eq. 4 makes clear the basic concepts for 
integrating the local weak form in Eq. 6. The MLPG formula- 
tion enables us to use different interpolations for the trial and 
the test functions. Furthermore, the sizes and shapes of the 
sub-domains of trial and test functions, respectively, do not 
need to be the same in the MLPG formulation. Therefore, the 
MLPG method can include all other meshless methods as spe- 
cial cases. In the present method, we use the same function 
space for the trial and the test functions as a special case. Note 
that the value of the trial function at each point x inside ft{, 
is influenced by a set of values of the function at an arbitrary 
number of nodes in the vicinity of each x, in a non-element, 
diffuse interpolation of the moving least square (MLS) type. 
Thus, Eq. 6 leads, for each ft{, to the I th equation in the system 


with 

— f E ( v for plane stress 

^ ~ | an V_ 1 *Pv for plane strain 

In the above equations, is the normal vector at the 

boundary, and E and v are the Young modulus and Poisson’s 
ratio, respectively. The local symmetric weak form in Eq. 9 
makes the "stiffness" entries, Kjj (which is the stiffness ma- 
trix in multi-dimensional space), in the row corresponding to 
the node /, and to the nodes J, depending only on the non-zero 
values of the integrands in the weak form, over the intersection 
of ft* and Then, the global equation can be written as 

y K,yU J = fy (10) 

J= l 

where 

K/v = [ (B[) T DB J dn 

JQI 
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«“(*,?) = p r (jt)a(jt) 



[/ 7 |(x),/? 2 (x)>--->Pm(x)] is a complete monomial basis 
of order m; and a(x) is a vector containing the coefficients 
a y -(x), j = 1,2,... which are functions of the space 
coordinates x = [*,}>, z] 7 . The commonly used bases in 2-D 
problems are the linear basis: 

p r (x) = [l.x.y] (12) 

or the quadratic basis: 

p T (x) = [l,x,y,x 2 1 xy,r] (13) 


Figure 3 : Conceptual explanation of the moving least square 
interpolation, in one-dimension. 

f ; = f v't</r+ [ v'mq 
J r', J n' 

Therefore, in the MLPG method, the usual assembly process 
is not required to form a global stiffness matrix. Theoretically, 
as long as the union of all local domains covers the global 
domain, the equilibrium equation and the boundary conditions 
will be satisfied in the global domain Q and its boundary T, 
respectively. 

3 Characteristics of moving least square (MLS) approxi- 
mation & numerical integration 


The coefficient vector a(x) is determined by minimizing a 
weighted discrete L 2 -norm, defined as: 


Y(x) = w,(x)[p 7 '(x,)a(x) - u‘ 


:/l2 


— [P- a(x) - u] r - W(x) [P- a(x) - u] 


(14) 


where w/(x) is a weight function defined in a sub-domain 
with the node / as its center; and with w/(x) > 0 for all x in 
the support of w/(x) and wy(x) = 0 at the boundary of x/ 
denotes the value of x at node /, and the matrices P and W are 
defined as 


r 


p = 


P' (Xl) 
P T (*2) 

P r (*w) 


(15) 


In this section, we review the characteristics of the moving 
least square (MLS) approximation, including some difficulties 
in the numerical integration of the weak form to evaluate the 
stiffness matrix, which was well explained in Atluri, Cho, and 
Kim (1999a) and Atluri, Kim, and Cho (1999b). The MLS 
method is generally considered to be one of the schemes to in- 
terpolate random data with a reasonable accuracy. Nayroles, 
Touzot, and Villon (1992) were the first to use the MLS in- 
terpolation in a Galerkin formulation, which they called the 
"diffuse element method". The MLS approximation always 
preserves completeness up to the order of the basis, and rea- 
sonably interpolates the irregularly distributed nodal informa- 
tion. However, the nodal shape functions that arise from the 
MLS approximation have a very complex nature. They are not 
only rational, but they are also of different types across the 
boundaries of sub-domains. This complexity results in diffi- 
culties with the numerical integration of the weak form in the 
MLPG method. 

We consider the approximation of a function u(x) in a local 
region centered at x in a domain Q as shown in Fig. 3. 

The moving least-square approximation starts from a local ap- 
proximation in the neighborhood of x, such as 


w- 


W|(x) 

0 


and 


0 


w w (x) 


J NxN 


~T r*i *2 
U 1 - [u\u-, 


,« w ] 


(16) 


(17) 


Here it should be noted that u 1 ,/ = 1,2, ■ • • , N are the fictitious 
nodal values. In fact, only those neighboring nodes 7, whose 
sub-domains intersect with the sub-domain Q l 5 of node /, 
have an influence in constructing the shape function for node 
/. For convenience, x in the above relations is replaced by x, 
because a local point x can be extended to all points in whole 
domain. This is the basic concept of the "moving” procedure, 
and we can finally obtain a global approximation. The concep- 
tual explanation for the MLS is given in Fig. 3. 

The stationary condition of F(x) with respect to the coeffi- 
cients a(x) leads to the following relation: 


A(x)a(x) = B(x)u 


(18) 


where the matrices A(x) and B(x) are given by 


u local {\,x) = p r (x)a(x) Vx € B(x) (1 1) 

where, B(x) is a sphere centered at x, p 7 (x) = 


A(x) = P T WP = jf w,(x)p(x i )p r (x,) 


/=! 


(19) 
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B(x) = P r W = [vvi(x)p(xi),W 2 (x)p(x 2 ),..-, 

w N (\) p(x w )] (20) 

The global approximation ^(x) can then be expressed as 
N 

u h (x) = '£«>' (x)u' (21) 

/= l 

where the nodal shape functions are given by 

<l> r (x) = p r (x) A -1 (x)B(x) (22) 

In the traditional Galerkin FEM, the ‘nodal shape functions’ 
have a value of unity at the respective node, and an approxima- 
tion of the type of Eq. 21 would involve the directly the ‘nodal 
value 7 of the field variable. However, in the present MLS ap- 
proximation, u l are fictitious, and are not exactly equal to the 
nodal values of the field variable (see Fig. 3). Inspite of this, it 
is instructive to call ^(x) in Eq. 21 ‘a nodal shape function’. 
The MLS interpolation is well defined only when the matrix 
A is non-singular. A necessary condition for a well-defined 
MLS interpolation is that at least m weight functions are non- 
zero (i.e. N > m) for each sample point x £ Q. The partial 
derivative of tf) 7 (x) can be obtained as follows: 

m 

<*= Sby^A-'B^ + p^A-'B.i + A^'B),/] (23) 

in which A"^. 1 = (A“ [ )^. represents the derivative of the inverse 
of A with respect to x, which is given by 

A-'zr-A-'A.i.A-' (24) 

and the index following a comma indicates a spatial deriva- 

tive. Considering that <|) / (x) = 0 whenever vv/(x) = 0, the sup- 
port sizes for the nodal shape function and the weight function 
have the same value. The nodal shape functions obtained by 
the MLS interpolation with mih order basis can reproduce any 
mth order polynomials g(x) exactly [Belytschko, Krongauz, 
Organ, Fleming, and Krysl (1996)], i.e., 

£= ♦'(*)*(*/)=*(*) (25) 

7=1 

Eq. 25 indicates that the nodal shape function is complete up 
to the order of the basis. In order to guarantee the convergence 
of the weak formulation with successive increase in the num- 
ber of nodes, the shape functions have to be complete. The 
smoothness of the nodal shape function ^(x) is determined 
by those of the basis and of the weight function. The choice 
of the weight function is more or less arbitrary as long as the 
weight function is positive and continuous. 

We can obtain an explicit form for the nodal shape functions, 
with a linear basis, in a two-dimensional problem, in order to 



Figure 4 : The nodal shape function <j) 7 (x) , from the MLS in- 
terpolation, which is nonzero in and which has a different 
form in each small intersection, as divided by several circular 
sub-domains. 


better understand the characteristics of the nodal shape func- 
tions. The nodal shape functions can be written as 

= '^^[c 0 {x,y),c [ {x,y),c 2 [x,y)} j * j (26) 

where the coefficients co(x,y), c \ (*>>’)> and b(x,y) are 

given in Atluri, Kim, and Cho (1999). In general, co(*,y), 
ci (x,y), C 2 (x,y), and b(x } y) are not of the single type of func- 
tions, because the sub-domains related to the sub-domain Q ! s 
make complex intersections as shown in Fig. 4. For exam- 
ple, the weight functions in the expressions for the coefficients 
co(x,y), ci(x,y), c' 2 (x,y), and b{x,y) are different at Xi and x? 
in Fig. 4. In other words, the nodal shape function (^(x con- 
sists of a different form of a rational function in each small 
intersection region, as indicated in Fig. 4. 

Although the smoothness of (^(x) can be achieved if a suffi- 
cient order of spline function is used as a weight function, the 
shape function in the sub-domain consists of several types 
of rational functions. It seems to be difficult to integrate these 
kinds of complex functions in the sub-domain, by using a sim- 
ple Gaussian quadrature rule, and this causes a difficulty in 
the numerical integration of weak forms. As a result, an ac- 
curate integration of the shape functions for the construction 
of the stiffness matrix is not as trivial as for the finite element 
method. Due to different characteristics of functions in differ- 
ent sub-regions, we cannot expect accurate numerical integra- 
tion even though many numbers of integration points are intro- 
duced. Two typical integration methods in meshless methods 
are illustrated in Fig 5. Although integration over sub-domain 
Q l s (Fig. 5b) is more natural than the integration using a back- 
ground cell (Fig. 5a), the difficulty of numerical integration 
may not be avoided in a simple way. 
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Figure 5 : Integration methods in meshless methods: (a) using 
background cell and (b) integration over sub-domain Q. ! s for 
node I. 


To surmount this difficulty, in the present version of the MLPG 
approach, polygonal sub-domains are employed so as to lead 
to a single type of the MLS shape function all over the inte- 
gration domain. The detailed explanation on how to construct 
the MLS shape functions, in the polygonal sub-domain Q* s , is 
given in the following section. 

4 Error control in the MLPG method 

This section first introduces the concept of primary and sec- 
ondary nodes, and describes the way to construct polygonal 
sub-domains and the corresponding weight functions, for 
the present version of the MLPG method. Later, the advan- 
tages of the present approach in the treatment of essential 
boundary conditions, and non-convex boundaries, are pointed 
out. Next, the simplified procedures for the numerical evalua- 
tion of the stiffness matrix is described. 

4.1 Polygonal sub-domains and weight functions for pri- 
mary and secondary nodes 

As explained in the previous section, the difficulty of numer- 
ical integration of the weak form is due to the complexity of 



(a) 



Figure 6 : Schematic representations for sub-domains with 
randomly distributed nodes: (a) for the MLPG method with 
circular sub-domains and (b) for the MLPG method with 
polygonal sub-domains; in (b), the sub-domain for the sec- 
ondary nodes is taken to be the same as that for the nearest 
primary node in the present method (The solid circles are pri- 
mary nodes; and the open circles are secondary nodes. 


the MLS shape functions in an integration domain. To evalu- 
ate the stiffness matrix with accuracy, it is necessary to make 
the shape functions simple in an integration domain. In our 
approach, we introduce the concept of primary and secondary 
nodes in the MLPG method, in order to prevent the crossing 
of the boundaries of sub-domains in an integration domain, 
such that the polygonal sub-domains defined by a mesh from 
only the primary nodes also become the sub-domains for the 
secondary nodes. As a result, no additional mesh is required 
for the secondary nodes, which makes it possible to extend the 
original MLPG concept to be a useful tool for error control 
and adaptive calculation. The shape functions in a sub-domain 
in the present version of the MLPG method may have a sim- 
pler form than those in other meshless methods, because of 
the alignment of the boundaries of sub-domains. Fig. 6 show's 
two types of sub-domains in the MLPG method: circular sub- 
domains and polygonal sub-domains that are used by both the 
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Figure 7 : Arbitrary secondary-node placement in the MLPG; 
(a) primary nodes, and (b) primary and secondary nodes in the 
domain. Primary nodes are solid circles, and secondary nodes 
are open circles. 


primary and the secondary nodes. Since the types of the shape 
functions may change across the boundaries of sub-domains, 
the first method shown in Fig. 6a may lead to great difficul- 
ties in integrating the weak form accurately, for randomly dis- 
tributed nodes, due to complex intersections of sub-domains. 
However, the present method in Fig. 6b uses simple polygo- 
nal intersections of sub-domains, in order to avoid this diffi- 
culty. As a result, the shape functions in a polygonal intersec- 
tion have a single type of a rational function, because of the 
alignment of the boundaries of sub-domains. The advantage 
of the second method, in Fig. 6b, is that the secondary nodes 
can be placed at any arbitrary locations, where it is needed to 
improve the accuracy of the solution, as shown in Fig. 7. This 
second method, which is attractive for adaptive error-control 
algorithm, is pursued in the present paper. Fig, 8 illustrates the 
comparison between finite element mesh/background mesh 
that is commonly used in other meshless methods such as 
EFG, RKPM and h-p cloud methods; and a background mesh 
of the primary nodes only as in the present method. It is impor- 
tant to note that in the present MLPG method, the secondary 
nodes do not involve an additional mesh. In addition, the sec- 
ondary nodes can be added and/or moved without changing a 
coarse background mesh constructed from the primary nodes. 
Fig. 8a shows a complex and a distorted feature for randomly 
distributed nodes in the FEM/EFG/RKPM/h-p cloud methods, 
and it is to be noted that a new mesh is required to change 
nodal information. 

To meet the staled goal, it is important to develop simple 
forms of the shape functions for the primary and the sec- 
ondary nodes. The weight functions in the MLS approxima- 
tion, which ultimately govern the shape functions (see Eqs. 
19, 20, and 22), should be defined appropriately to generate 
shape functions with simple forms. In the present approach, 
the secondary nodes can be placed at arbitrary locations in the 
domain Q. and on the global boundary T, in order to improve 
the deformations in each region, respectively. First, the sec- 
ondary nodes which may be placed randomly, in the domain 



Figure 8 : Comparison between the FEM mesh/background 
mesh, and the MLPG: (a) the mesh from the usual Finite Ele- 
ment Method as well as the background mesh in EFG, RKPM, 
and h-p cloud method, when primary and secondary nodes are 
used, and (b) the mesh from primary nodes in the MLPG ap- 
proach, whereas the secondary nodes in the present method do 
not involve an additional mesh. 



Figure 9 : A primary-node anchor for a secondary node is the 
nearest primary node. /, 7, K and L indicate the primary nodes, 
and i, j, k and / indicate the set of secondary nodes. 
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centered at the secondary node / 

Figure 10 : The definition of sub-domains for the primary 
node I and the secondary node i in the domain £1 r, is the ra- 
dius of support for spline function to construct weight function 
for the secondary node / in the domain. 


Q, take the sub-domains to be the supports of the nearest pri- 
mary nodes, as shown in Fig. 6. The primary nodes on the 
global boundary F are excluded from the set of candidate pri- 
mary nodes that act as anchors for the secondary nodes in the 
domain Q, so that the essential boundary conditions can be im- 
posed directly on the primary nodes on the global boundary F. 
Fig. 9 indicates how to choose the primary nodes that would 
be anchors for the secondary nodes in the domain, in a simple 
way. Therefore, it is very easy to find the sub-domains for the 
secondary nodes in the domain, with an initial coarse back- 
ground mesh, constructed only from the primary nodes. We 
denote by I the primary node and by i the set of the secondary 
nodes, related to the primary node /, as shown in Fig. 10. 
There are several methods to construct the weight functions on 
polygonal sub-domains. In the present paper, we use the linear 
finite element shape functions as weight functions in the MLS 
approximation for primary nodes, i.e., 

w,(\) = N'(x) on Q' (27) 



weight function for weight function for 

primary node / secondary node i 

(a) 


shape function for shape function for 

primary node / secondary node / 

(b) 

Figure 11 : Weight and shape functions: (a) weight functions 
and (b) shape functions for the primary node / and for the sec- 
ondary node /, respectively in the domain. The weight func- 
tion for the primary node is the finite element shape function, 
and the weight function for the secondary node has a skewed 
form which is obtained by multiplying the finite element shape 
function centered at /, and the 4 th order spline function cen- 
tered at i. 

at the secondary node /, as indicated in Fig. 10. Within this cir- 
cular support centered at the secondary node i, the maximum 
value of the spline function is at the position of the secondary 
node /, and the size of the circular support is set by the maxi- 
mum distance from the secondary node i to the primary nodes 
in the sub-domain Q*. We take the weight function for the 
secondary node i to be 




where A^(x) is the linear finite element shape function for the 
primary node I. Note that the weight function inside a cell 
all over the sub-domain Q f s has only one type of a simple 
form. We now discuss the construction of the weight func- 
tions in the MLS approximation for the secondary nodes. The 
weight functions for the secondary nodes should be zero at 
the boundaries of polygonal sub-domains as defined for the 
original primary nodes, which is the fundamental requirement 
in the MLPG formulation. For this purpose, we choose the 
weight functions in the MLS approximation, for the secondary 
nodes, to be the product of the linear finite element shape func- 
tion N ! (x) centered at a primary node / as in Eq. 27, and the 
4 th order spline function w,(x) with a circular support centered 


Wj(x) = N 1 (x)\Vj{x) on = Q.‘ s (28) 

with 

1 -6(f) 2 + 8(^) 3 -3(f)\ 0 <d,<n 

0, d, > r, 

(29) 

where I and i are the primary and the secondary nodes in the 
domain Q, respectively, = jx- x,-| is the distance from the 
secondary node x,-, and the radius r t is the size of the support 
chosen in such a way that the circular support for w,(x) covers 
the polygonal sub-domain Q[ — Note that the values of 
weight functions in Eq. 28 for the secondary nodes are zero 
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support for the spline function vvy(je) 
centered at the secondary node i 


Figure 12 : The definition of the sub-domains for the primary 
node I and the secondary node i on the global boundary V . 
r, is the radius of support for spline function to construct the 
weight function for the secondary node L 


at the boundary of the sub-domain £2^ = £)*. The weight and 
the shape functions for the primary and the secondary nodes 
in the domain are plotted in Fig. 11. In these figures, the 
weight function for the secondary node i has a skewed form, 
depending on the location of x,-. The construction of the shape 
functions and derivatives for the secondary nodes is identical 
to the procedure used for the primary nodes. The weight func- 
tions for secondary nodes in the domain satisfy the following 
conditions: 

w/(x) > 0 on the sub-domain £}' = 

vv,(x) = 0 at the boundary of the sub-domain £l l s = 

H\(x) has a simple form of a continuous function inside a cell 
all over the sub-domain Q. l s = £^. 

The deformations on the global boundary cannot be improved, 
by adding secondary nodes only in the domain, because the 
primary nodes on the global boundary are excluded as can- 
didate primary nodes for the secondary nodes in the domain. 
This limitation may spoil the quality of numerical solutions, 
due to the limitation of only a linear deformation between the 
primary nodes on the global boundary. However, improved 
deformations on the global boundary are important, in order 
to obtain the stresses or the strains correctly on the global 
boundary. Also, the convergence of energy may not be at- 
tained, without an improvement of the deformations on the 
global boundary, even though the convergence of displace- 
ments may be obtained using only the secondary nodes in the 
domain. Hence, in order to improve the deformations on the 
global boundary, some secondary nodes are placed also on 
the boundary-segments connecting the primary nodes on the 
global boundary. The sub-domains for the secondary nodes on 
the global boundary are the cells connecting primary nodes on 



weight function for weight function for 

secondary node / secondary node j 


(a) 



shape function for shape function for 

secondary node i secondary node j 

(b) 


Figure 13 : (a) weight functions and (b) shape functions for 
secondary nodes i and j on the global boundary T. 

the global boundary, as shown in Fig. 12. The weight func- 
tions for the secondary nodes on the global boundary can be 
constructed as 

vv,(x) = jY m (x)w/(x) in £2^ (30) 

where w,(x) is the 4 ,h order spline function in Eq. 29, where 
r, is the maximum distance from the secondary node / to the 
primary nodes in the sub-domain Q^, and /V m (x) in the present 
method is taken by 

Ar(x) = ^M = I(i-^ 2 )(i-n) (3i) 

In the above equation, the E, and rj are the coordinates defined 
on the master domain, as in the finite element method. The 
weight and the shape functions for the secondary nodes on the 
global boundary are plotted in Fig. 13. In these figures, the 
secondary nodes / and j on the global boundary between two 
primary nodes have the same sub-domain Q l s and Q J S . Conse- 
quently, the linear deformations on the global boundary can be 
improved by adding the secondary node on the global bound- 
ary. 

To alleviate the difficulty in the numerical integration of the 
weak form, it is important to preserve a single type of a contin- 
uous function in an integration domain. Since the MLS shape 
functions, which are the rational functions as in Eq. 26, are 
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• primary node 
n secondary node 
on boundary 


Figure 14 : Schematic representation of the shape functions 
for the secondary nodes / and which are located between the 
primary node I and J on the global boundary F, 



Figure IS : Sub-domains for a domain with non-convex 
boundaries (a) for other meshless methods with circular sub- 
domains, and (b) for the present MLPG method. 


expressed in terms of the weight functions and the basis in the 
MLS approximation, the integrand in the weak form within an 
integration domain consists of a continuous rational function 
in the presently described formulation. Again, the emphasis is 
placed on the fact that the shape functions for the primary and 
the secondary nodes have a single form of a rational function 
all over the integration domain, because there is no crossing 
of the boundaries of the sub-domains in the domain for inte- 
gration. Although the MLS shape functions in a sub-domain 
become more complex as the number of the secondary nodes 
increases, a single form of a rational function may be much 
easier to integrate numerically, than the shape functions with 
different types of functions as in other meshless methods. 

In summary', the secondary nodes can be added at arbitrary 
positions in the domain and on the global boundary, after con- 
structing the initial polygonal sub-domains from the primary 
nodes, and errors in numerical results can be controlled by 
adding and moving the secondary nodes. Thus, to start with, 


C1 J S =Q 


i 

s 



Figure 16 : Numerical integration inside each integration cell 
in the sub-domain Q ! s = £3' . 


only a few primary nodes may be used, and later, a random 
pattern of the secondary nodes may be introduced, in an adap- 
tive fashion to control the error of the numerical solution. 


4.2 Treatments of boundary conditions and non-convex 
boundaries 

One major difficulty in the meshless methods is considered 
to be the imposition of the essential boundary conditions, be- 
cause, in general, the approximation functions do not satisfy 
the Kronecker-delta condition ty l (xj) = at the boundary. 
Most meshless methods have used Lagrange multipliers or 
penalty methods to impose the essential boundary conditions 
[Zhu and Atluri (1998)]. In some cases, meshless interpo- 
lations and FEM shape functions have been combined [Be- 
lytschko, Organ, and Krongauz (1995)], leading to a complex 
interface element in the regions of intersection of FEM and 
meshless shape functions. However, in the present method 
the Kronecker-delta condition is satisfied at the primary nodes 
on the global boundary, because the primary nodes on the 
global boundary are not candidates as anchors for the sub- 
domains for the secondary nodes in the domain, and because 
the weight functions, except those for the primary nodes on the 
global boundary, are zero at this point. Therefore, the essen- 
tial boundary conditions can be imposed exactly at the primary 
nodes on the global boundary. On the contrary, the Kronecker- 
delta condition may not be satisfied at the secondary nodes 
on the global boundary, as a result of non-zero values for the 
weight functions at the secondary nodes. Since the deforma- 
tions on the global boundary depend only on the primary and 
the secondary nodes on the boundary-segment between the pri- 
mary nodes on the global boundary, the fictitious nodal values 
for the secondary nodes on the global boundary can be evalu- 
ated in an easy way. Fig. 14 illustrates two secondary nodes 
on the boundary-segment between two primary nodes on the 
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Figure 17 : Geometric description for numerical experiments: (a) a cantilever beam, (b) a plate with a hole, (c) concentrated 
load on a semi-infinite plate, and (d) a center cracked plate in tension. 


global boundary, and the following equation can be written as 


4>'(x/) <t> 7 jx/) <t> y (x/) ' 


u> ] 


u[ \ 

0'(x/) <t>'(x/) 4> ; (x/) <t > J (X/) 

J 

u‘ 


Ti l 1 

<t»'(x/) 4>'(x/) <!> 7 (x/) <t> y (x/) 

\ 

I uJ I 

► — ^ 

W [ 

<t»'(x,) <|>''(x/) <j)^(x,) (t> y (x/) 


l J 




(32) 

where u is the fictitious nodal displacement in MLS approx- 
imation, and u is the prescribed displacement on the global 
boundary. As explained before, the following conditions are 
satisfied at the primary nodes: 

u 1 = u' and ii J = r/ (33) 

We can evaluate the fictitious nodal values u l and u J at the 
secondary nodes on the global boundary from Eq. 32, and 
the essential boundary conditions can be imposed directly in 
the computation. Consequently, we can impose the essential 
boundary conditions exactly at the primary and the secondary 
nodes on the global boundary. 

The other meshless methods with circular sub-domains may 
lead to a difficulty near boundaries when the domain is 


not strictly convex [Organ, Fleming, Terry, and Belytschko 
( 1 996)]. In Fig. 1 5a, the deformation at the node / on the upper 
part of edge directly affects the deformations on the lower part, 
which is invisible to an observer at the node /. Some special 
treatments, as introduced by Organ, Fleming, Terry, and Be- 
lytschko (1996) to solve the problems with non-convex bound- 
aries, are required in other meshless methods. In the present 
method, however, there is no this difficulty with a non-convex 
boundary, because the sub-domains for the primary and the 
secondary nodes are simple polygonal-local domains as shown 
in Fig. 15b. Therefore, we can easily deal with problems with 
non-convex boundaries using the present method. 

4.3 Numerical integration of weak forms 

To evaluate the stiffness matrix from the weak form, it is nec- 
essary to use a numerical quadrature since analytical integra- 
tion is all but impossible in general. The numerical integra- 
tion of the stiffness matrix usually plays an important role 
in the convergence of numerical solutions in meshless meth- 
ods. Fig. 5 shows that the schematic features of two integra- 
tion methods in meshless methods. The first method, using 
a background mesh, has been used in most meshless meth- 
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Figure 18 : Comparison of the MLPG and the exact solutions, for the problem of a cantilever beam. Displacements and 
distribution of stress G\\ (a) without secondary nodes, (b) with 50 secondary nodes, (c) with 100 secondary nodes, and (d) with 
200 secondary nodes in the domain. 


ods including the EFG, RKPM, and hp-clouds. Dolbow and 
Belytschko (1999) have already indicated that the integration 
using the background mesh is not adequate, to accurately in- 
tegrate the terms in the stiffness matrix, when irregularly dis- 
tributed nodes are used. They presented a method to reduce er- 
rors in numerical integration, by making the integration cell to 
be aligned with the boundaries of sub-domains. As explained 
before, the present method, however, leads to a simple type of 
a rational function inside a cell in a sub-domain. Therefore, 
a cell (a quadrilateral for example) in a sub-domain {Q ! s ~ 
for the MLPG method) is taken as an integration domain in 
the present MLPG formulation, to obtain accurate numerical 
integration for the stiffness matrix. The schematic of the inte- 
gration method is presented in Fig. 16. 

Gaussian quadrature is commonly employed to numerically 
evaluate the integrals in the weak forms. The Gaussian quadra- 
ture can exactly integrate the polynomials of order 2/? — 1 in a 
spatial direction, where n is the number of integration points in 
that spatial coordinate. Since inside an integration domain the 


shape functions are rational functions in the present method, 
the Gaussian quadrature may not be adequate to evaluate inte- 
grals in the present method properly. However, the accuracy of 
numerical integration may be controlled by the number of in- 
tegration points, taking a proper level of polynomial as an ap- 
proximation for these rational functions. In general, the more 
secondary nodes are involved in the sub-domain, the more in- 
tegration points are required. To find an efficient integration 
rule for rational functions is still an open question, to improve 
the performance of this method while using only a small num- 
ber of integration points. 

5 Numerical experiments 

Several problems in two-dimensional linear elasticity are 
solved to illustrate the effectiveness of the present method. 
The numerical results of the MLPG method as applied to 
problems in two-dimensional elasto-stalics, specifically a can- 
tilever beam, a plate with a hole (circular and elliptic holes), 
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concentrated load on a semi-infinite plate, and a center cracked 
plate in tension as shown in Fig. 17, are now discussed. The 
Young modulus and the Poisson’s ratio are E = l.Ox 10 10 and 
v = 0.25, respectively. We use the displacement and energy 
norms defined as 


IMI=(/ u r ucK2)5 
II e IN (x f e r DeJQ)i 

2 Jn 

The relative errors for || u || and || e |j are defined as 
_ || u"""' - u" acl || 

"" “ II u«-' || 

jj £ num _ £ exact jj 

rt= ifi=nr 


(34) 

(35) 


(36) 

(37) 


The linear basis in the MLS approximation is used in the nu- 
merical examples. 


5.1 Cantilever beam 


We first consider a cantilever beam problem shown in Fig. 17a. 
The exact solution for this problem is given in Timoshenko and 
Goodier (1970) as 


Ml - j)[M^r-x)+(2 + v)y{y-D)] (38) 

" 2= db [ ^ (3L ~* ) + 3v(L-;c)(y_ f ) 2 + ^J^ 2 *) (39) 


where 


The stresses corresponding to the above are 


P D 

On ~--(L-x){y- -) 

022 — 0 

012 = -jj{y-D) 


(40) 

(41) 

(42) 


We use regularly distributed primary' nodes for a model with 
D — 4.0 and L = 8.0 to examine the effects of the secondary 
nodes. The essential and traction boundary conditions are ap- 
plied at the left and right sides of the beam, respectively, and 
P in Eqs. 38-42 is (— ) 1 .0 x 10 8 . The solution without the sec- 
ondary nodes is exactly the same as that of the finite element 
method, as shown in Fig. 18a. We distribute 50, 100, and 
200 secondary nodes randomly in the domain through gener- 
ating random numbers. To improve the deformations on the 
global boundary, the secondary nodes are placed between the 
primary nodes on the global boundary, as shown in Fig. 18. 
We use 3 x 3, 5 x 5, and 8x8 integration points in an inte- 
gration domain for 50, 100, and 200 secondary nodes in the 



Figure 19 : Convergence for the displacement and the en- 
ergy norms in a cantilever beam problem, with irregularly dis- 
tributed secondary nodes in the domain. 


domain, respectively. It is required to increase the number of 
integration points as the secondary nodes are added. Com- 
parison between the MLPG and the exact solutions arc given 
in Fig. 18. In these results, the relative displacement error 
changes from 9.26% to 0.38%, and the relative energy error 
decreases from 31.37% to 10.74% by adding 200 secondary 
nodes. In Fig. 19, the convergences in the relative displace- 
ment and energy norms for the MLPG are presented with a 
representative nodal density y/N~ 5 , where N s is the number of 
secondary nodes in the domain. The present method gives bet- 
ter results than the solutions without the secondary nodes. Par- 
ticularly, the relative displacement error decreases rapidly even 
though the secondary nodes are added randomly. 


5.2 Infinite plate with a hole 

We consider an infinite plate with a circular hole of radius a 
(b — a in Fig. 1 7b). The plate is subjected to a uniform tension, 
ao = 1 .0 x 10 9 , in the x-direction, at infinity as shown in Fig. 
17b, The exact solutions for stresses are 

cfi 3 3a 4 

an = Oo[l - ^ (-cos 20 + cos 40) + ^“4 cos 40] (43) 

(p- 1 3 a^ 

0 12 = Oo [-^2 (- sin20 + sin40) + — sin40] (44) 

c p 1 3 a^ 

a 22 = o 0 [ — T (- cos 20 -cos 40) - — T cos40] (45) 

n 2 2r a 

where (r, 0) are the polar coordinates, and 0 is measured from 
the positive x-axis in a counterclockwise direction. The cor- 
responding displacements, in the plane stress case, are given 
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Figure 20 : Comparison of the MLPG and the exact solutions, for the problem of a plate with a circular hole. Displacements 
and distribution of stress G\ \ (a) without secondary nodes, (b) with 50 secondary nodes, (c) with 100 secondary nodes, and (d) 
with 200 secondary nodes in the domain. 
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quadrant of the plate is modeled as shown in Fig. 17b. Sym- 
metry conditions are imposed on the left and bottom edges, i.e. 

~ 0, tj = 0 is prescribed on the left edge and «2 = 0, t\ =0 
on the bottom edge, and the inner boundary at a — 1 0 is trac- 
tion free. The traction boundary conditions, as given by the 
exact solutions, are imposed on the outer boundary at r — 4.0. 
Secondary nodes are distributed randomly in the domain, as 
shown in Fig. 20. 4 x 4, 6 x 6, and 9x9 integration points in 
an integration domain are used for models with 50, 100, and 
200 secondary nodes in the domain, respectively. One sec- 
ondary node, on each boundary-segment between the primary 
nodes on the global boundary, is added to improve the defor- 
mations on the global boundary. The comparison between the 
MLPG and the exact solutions are shown in Fig. 20. 


Due to symmetry, only a part, 0 < r < 4, of the upper right 
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Figure 21 : Convergence for the displacement and the energy 
norms in the problem of a plate with a circular hole, with ir- 
regularly distributed secondary nodes in the domain. 


The relative displacement error decreases from 5.93% to 
0.48%, and relative energy error decreases from 11.51% to 
4.69% by adding 200 secondary nodes in the domain. The 
convergence with the number of nodes is shown in Fig. 21 
with a representative nodal density y/N~ s for the irregularly dis- 
tributed secondary nodes in the domain. 

Agreement between the MLPG and the exact solutions is ex- 
cellent in this example. As a more severe case, we consider a 
plate with an elliptic hole for aspect ratios of 4 : 1 , and 8 : 1 , re- 
spectively. The stress concentration factors for these elliptical- 
hole problems in an infinite plate under uniform traction are 
9.0 and 17.0, respectively. The secondary nodes in the domain 
and on the global boundary are added near the elliptical hole 
where the major-axis is of the ellipse intersects the hole sur- 
face. The models with randomly distributed (25, 50, 100 and 
200) secondary nodes in the domain, use 4x4, 6x6, 9x9 and 
11x11 integration points in an integration cell, respectively. 
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(a) 


sii 

9.0000E+09 
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7.7143E+09 
7.071 4E+09 
6.4266E+09 
5.7857E + 09 
5.1429E+09 
4.5000E+09 
3.8571 E+09 
3.2143E+09 
2.571 4E +09 
1 .9286E+09 
1.2857E+09 
6.4286E+08 
OOOOOE 1 00 


plate with an el&ptic hole (b/a = 8.0) 
number of secondary nodes in domain = 25 
* primary node 
c secondary node in domain 
■ secondary node on boundary 



maximum stress = 1 .0681 E + 1 0 
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plate with an eBiptic hole (b/a = 4.0) 
number ol secondary nodes in domain = 50 
• primary node 
0 secondary node in domain 
. secondary node on boundary 





S11 

9.0000E+09 
8.3571E+09 
7.7143E+09 
7.071 4E+09 
6.4286E+09 
5.7857E+09 
5 1429E+09 
4. 5000 E +09 
3.8571 Ef 09 
3.2143E+09 
2.571 4E+09 
1 .9286E+09 
1.2857E+09 
6.4286E+08 
O.OOOOE+OO 


maximum stress = 8.8231 E+09 
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O.OOOOE+OO 


maximum stress = 1 .386 1 E + 1 0 


(C) 


plate with an elliptic hole (b/a = 4.0) 
number of secondary nodes in domain - 1 00 

* primary node 

o secondary node in domain 

• secondary node on boundary 


o* 

S? 80 

V ' 


•oo 
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I 9.0000E+09 
8.3571 E+09 
7.7143E+09 
7.071 4E +09 
6.4286E+09 
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5.1429E+09 
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3.2143E+09 
2.5714Et09 
1 9286E+09 
1 .2857E+09 
6.4266E+08 
O.OOOOE+OO 


maximum stress = 9.081 5E+09 


plate with an eBiptic hole (b/a = 8.0) 
number of secondary nodes in domain = 100 
. primary node 
e secondary node in domain 
■ secondary node on boundary 
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• o * C o* 

A ^ o co # ©o 
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1 7000E » 10 
1.5786E+10 
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6.071 4E » 09 
4.8571 Ef 09 
3 6429E.09 

2 4286E » 09 
1.2143E+09 
O.OOOOE+OO 


maximum stress - 1 .4645E+10 


(d) 


plate with an elSplic hole (b/a = 4.0) 
number of secondary nodes in domain = 200 
• primary node 
<> secondary node in domain 
■ secondary node on boundary 
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7.7143E»09 
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4.5000E+09 
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3.2143E+09 
2.571 4E+09 
1.9286E + 09 
1 2857E+09 
6.4266E + 08 
O.OOOOE+OO 


maximum stress - 9.0621 E+09 


plate with an elliptic hole (b/a = 8.0) 
number of secondary nodes = 200 
• primary node 
o secondary node in domain 
. secondary node on boundary 



Sll 

1 7000E+10 
1 5786E+ 10 
1 .457 1 E ♦ 1 0 
1.3357E+10 
1.2143E+10 
1 0929E+ 10 
9.7143E+09 
8.5000E 1 09 
7.2857E 1 09 
6.0714E+09 
4.8571 E. 09 
3 6429E+09 
2.4286E+09 
1.2143E+09 
O.OOOOE+OO 


maximum stress = 1.6437E+10 


(e) 


Figure 22 : The MLPG solutions for the problem of a plate with an elliptical hole: (b/a — 4.0 and b/a = 8.0). Displacements 
and distribution of stress C\\ (a) without secondary nodes, (b) with 25 secondary nodes, (c) with 50 secondary nodes, (d) with 
100 secondary nodes, and (e) with 200 secondary nodes in the domain. 
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Figure 23 : Convergence for the displacement and the energy 
norms in the problem of a plate with an elliptical hole, with 
irregularly distributed secondary nodes in the domain. 


Note that the sub-domains not influenced by the secondary 
nodes require only 2x2 integration points because of sim- 
ple shape functions in this region. Fig. 22 shows the stress 
distributions by adding 25, 50, 100 and 200 secondary nodes. 
The stress On near the tip of the major-axis of the elliptical 
hole is plotted in Fig. 23 with a representative nodal density 
y/N~ s , The results show that the stress concentration factors ap- 
proach the exact values as the number of the secondary nodes 
increases. 


5.3 Concentrated load on a semi-infinite plate 

Consider the concentrated load on an infinite plate as shown 
in Fig. 17c. The plate is assumed to be of unit width so that 
the concentrated load is equal to q. The exact solution for this 
problem [Saada (1974)] is given by 


O rr = 


2 q cos G 
71 r 


oee = o r e = o zz — 0 


e rr = 


t'ee = 


2 q cosG 
71 E r 

2 vq cos G 
tiE r 


(48) 

(49) 

(50) 

(51) 


ere - 0 

The displacements are 

2q ^ , d (1 — v)^ . 

u r — — ^cosG in — 0sm9 

nE r nE 

o(I-l-v) . ^ 2q . d (1— v)^ 

mq = — ^ sin 0 si n 0//7 — -GcosO 

nE nE r nE 


(52) 


(53) 

(54) 


in Fig. 17c. The concentrated load q is 1.0 x 10 9 in this ex- 
ample. Since the displacements are singular at the point of 
application of the concentrated load q, we use the model with 
a small circular hole a = 0.5. The symmetric boundary condi- 
tion is applied at the left side of the half model, and tractions 
are applied at the other boundaries from the exact solution. In 
numerical calculation, the distance d is set to be 4.0 to give 
the fixed boundary condition at the center. Convergence stud- 
ies are carried out using three different numbers of secondary 
nodes, namely 50, 100, and 200. 4x4, 6x 6, and 9x9 inte- 
gration points are used in an integration domain for these three 
models. The comparison between the MLPG and the exact so- 
lutions is plotted in Fig. 24, and the convergence is shown in 
Fig. 25. In this example, we obtain the similar trends as shown 
in the problems of a cantilever beam and a plate with a hole. 

5.4 Center cracked plate in tension 

Next, we consider a center cracked plate in tension. Due to 
symmetry, the right half as shown in Fig. 17d is modeled un- 
der plane stress condition. The size of model is h = w = 4 0, 
and the crack length is a = 2.0. The applied stress Go at the 
top and the bottom is 5.0 x 10 8 in this example. The symmetric 
condition is applied on the left side. Of primary importance in 
a crack problem is the determination of the parameters which 
characterize the singularity of the stress fields in the vicinity of 
a crack tip. The mode I stress intensity factor K}, as a charac- 
terizing parameter for the crack, is computed from 7-intcgral 
using domain integration [Nikishkov and Atluri (1987) and 
Anderson (1991)]. The size of the 7-integral domain is cho- 
sen as 2b\ x 2b2 = 2.0 x 2.0. The stress intensity factor Kj 
is evaluated by Ki = y/TE for plane stress, and the target so- 
lution for this problem is Ki/Kq = 1 .325 where Kq ~ Gq y/na 
[Tada, Paris, and Irwin (1977) and Wu and Carlsson (1991)]. 
In this numerical example, we use regularly distributed pri- 
mary nodes as shown in Fig. 26, and secondary nodes are 
added randomly near the crack tip. 4x4, 6x6 and 8x8 
integration points are used for models with 25, 50 and 100 
secondary nodes in the domain, respectively. However, only 
2x2 integration points are used in the sub-domains which are 
not influenced by the secondary nodes. The stress G 22 at the 
crack tip is higher as the number of the secondary nodes near 
the crack tip increases. By adding secondary nodes on the 
crack lines, the improved deformations near the crack tip are 
obtained as shown in Fig. 27. The deformations of the crack 
tines are very important to recover the stresses correctly near 
the crack tip. In Fig. 28, the errors in the stress intensity fac- 
tors evaluated from 7-integral are plotted against the number 
of the secondary nodes. The stress intensity factor approaches 
the target solution as the number of secondary nodes in the 
domain increase. 


where d is a distance from the point of the concentrated load 
q . We use the right half model with 4.0 x 4.0 size as shown 
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number of secondary nodes in domain = 0 


. primary node 



relative displacement error = 0.09264 
relative energy error = 0.2844 


number of secondary nodes in domain = 50 
• primary node 


° secondary node in domain 



relative displacement error = 0.0352 
relative energy error = 0.1 649 


(a) 


(b) 


number of secondary nodes in domain = 100 
. primary node 
* secondary node in domain 



number of secondary nodes in domain = 200 
. primary node 


° seconda ry node in doma in 



relative displacement error = 0.0209 relative displacement error = 0.01 59 

relative energy error = 0.1 347 relative energy error = 0.1 1 87 


(c) (d) 

Figure 24 : Comparison of the MLPG and the exact solutions for the problem of a concentrated load on a semi-infinite plate. 
Displacements and distribution of stress G 22 (a) without secondary nodes, (b) with 50 secondary nodes, (c) with 100 secondary 
nodes, and (d) with 200 secondary nodes in the domain. 


6 Concluding remarks 

A method, based on the meshless local Petrov-Galerkin 
(MLPG) concept for solving boundary value problems has 
been presented in this paper. The concept of primary and sec- 
ondary nodes was introduced, and the weight functions for 
primary and secondary nodes were defined to yield a single 
type of a rational function all over the integration domain. The 
approach presented in this paper alleviates a major difficulty 
in the numerical integration to evaluate weak forms in other 
meshless methods. Additional mesh for the secondary nodes 
is not required to improve the accuracy of solution. The essen- 
tial boundary conditions are easily taken care of as in FEM, 


and the non-convex boundaries can be treated without a spe- 
cial technique. 

A clear advantage of the present method is that the secondary 
nodes can be placed at any random locations, without the bur- 
densome task of constructing a new mesh to enrich the so- 
lution. The present MLPG method can control errors in nu- 
merical results by adding the secondary nodes, using only 
the sub-domains constructed from the primary nodes. There- 
fore, the present method can be a useful tool for error con- 
trol and adaptive calculation in the field of computational me- 
chanics. The numerical experiments show that the present 
method is very efficient for adaptive error control by placing 






Arbitrary placement of secondary nodes, and error control, in the Meshless Local Petrov-Gahrkin (MLPG) method 


29 



Figure 25 : Convergence for the displacement and the energy 
norms in the problem of a concentrated load on a semi-infinite 
plate, with irregularly distributed secondary nodes in the do- 
main. 


secondary nodes arbitrarily in the domain and on the global 
boundary. The present method can be easily implemented in 
three-dimensional problems with the same advantages. Fur- 
ther results based on the current approach will be presented in 
a series of forthcoming papers. 
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center cracked plate 

number of secondary nodes in domain = 0 



primary node 
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center cracked plate 

number of secondary nodes in domain = 25 



primary node 

secondary node in domain 
secondary node on boundary 



S22 

2.0000E+09 
1 .8571 E+09 
1.7143E+09 
1.5714E+09 
1 4286E+09 
1 .2857E+Q9 
1 .1429E+09 
1 .0000E+09 
8.571 4E+08 
7.1429E+08 
5.7143E+08 
4.2857E+08 
2.8571 E+08 
1 .4286E+08 
0.0000E+00 


Kl / KO = 1 .2647 (4.5 % error) 


(b) 


center cracked plate 

number of secondary nodes in domain = 50 


center cracked plate 

number of secondary nodes in domain = 100 


• primary node 
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secondary node in domain 
secondary node on boundary 
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Kl / KO = 1 .3107 (1 .08 % error) 


Kl / K0 = 1 .3330 (0.604 % error) 


(c) 


(d) 


Figure 26 : The MLPG solutions for the problem of a center cracked plate. Displacements and distribution of stress O 22 (a) 
without secondary nodes, (b) with 25 secondary nodes, (c) with 50 secondary nodes, and (d) with 100 secondary nodes in the 
domain. 
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Figure 27 : Deformation near crack tip, including secondary 
nodes on the crack line. 



sqrt (number of secondary nodes in domain) 

Figure 28 : Errors in the evaluation of stress intensity factor, 
versus the number of secondary nodes in the domain for a cen- 
ter cracked plate problem. 



